Transitions to Nematic states in homogeneous suspensions of high aspect ratio 

magnetic rods 



Gopinath A.|, Mahadevan L.|, and Armstrong R. C.,§ 

| Division of Engineering and Applied Sciences, 

Harvard University, Cambridge, MA 02138 

^Department of Chemical Engineering, MIT, Cambridge, MA 02139. 

(Dated: February 2, 2008) 

Abstract 

Isotropic-Nematic and Nematic-Nematic transitions from a homogeneous base state of a suspension of 
high aspect ratio, rod-like magnetic particles are studied for both Maier-Saupe and the Onsager excluded 
volume potentials. A combination of classical linear stability and asymptotic analyses provides insight into 
possible nematic states emanating from both the isotropic and nematic non-polarized equilibrium states. 
Local analytical results close to critical points in conjunction with global numerical results (Bhandar, 2002) 
yields a unified picture of the bifurcation diagram and provides a convenient base state to study effects of 
external orienting fields. 
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Recently, a kinetic theory based model for dispersions of acicular magnetic particles was 
developed 1 ' 2 using ideas grounded in classical models for liquid-crystalline polymers 3 . Effects of 
Brownian motion, anisotropic hydrodynamic drag, a steric force chosen to be of the Maier - Saupe 
form and a mean-field magnetic potential were included. Both continuum descriptions obtained 
via closure approximations and the diffusion equation were solved numerically for some parameter 
ranges 1,2 . The focus of this article is on obtaining a theoretical characterization of transitions to 
nematic states from a homogeneous base state of a suspension of slender high aspect ratio mag- 
netic particles. Combining local asymptotic and stability analysis near critical points with global 
numerical results, we obtain a physically convenient point of departure for investigations of ex- 
ternal aligning fields. Both the Maier-Saupe and the Onsager potentials are considered. Results 
for the Maier-Saupe potential are in excellent agreement with available numerical solutions of the 
equations and complement recent investigations on the classical Doi model 4 . 

The particles in the homogeneous dispersion are modeled as two point masses connected by 
a rigid massless rod of length L and diameter d with inherent magnetic dipoles, the magnetic 
moment being along the axis 1 ' 2 . We envisage a situation in which d and L are kept constant 
and the concentration of the rods can be varied. The orientation of the rod is specified by the 
unit vector u along the axis from one specified bead to another. In the mean-field approximation 
it suffices to consider one test particle in a sea of others. Denoting the orientation distribution 
function by f(u,t), one writes for the case of constant diffusivity in scaled form 6 

^=Stu-(Xuf + Wu(VEV + V M )). (1) 

Here 5ft u (-) is the rotation operator and the potentials are measured in units of fc^T. We define 
the average of a quantity, X(u), as (X(u)) = J X(u) /(u)du. The excluded volume intermolecular 
potential for a Maier-Saupe (MS) or Onsager (O) potential can then be written as 

V EV (u) = J (3 MS/0 (u, u') /(u', t) du', (2) 

where, /3ms(u, u = —^-ms( u ■ u ') 2 > ^ms being a phenomenological constant proportional to the 
concentration of rods, N and /3o(u, u') = 2NL 2 d |u x u'|. The total potential due to the mean 
magnetic field, Vm, can be written 2 

V M = -(3/2)£'<uu> : uu - A'u • (u) +A a + B (3) 

The first term reflects a net magnetic interaction potential due to average order 1 ' 2 , the second term 
is the mean field approximation to the dipole-dipole interaction between particles and A and B Q 
are constants independent of u. 
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Equations (l)-(3) do not involve any preferred direction for orientation of possible nematic 
states and so we choose to employ an expansion for /(u, t) in terms of spherical harmonic functions 
Yj m (u) = Y l m (9, <p). where u = (sin 9 sin 4>) e x + (sin 9 cos <p) e y + (cos 9) e z and e z is the axis from 
which 9 is measured. Since / is real valued, we can write 

oo +1 

/(«.*) = E E b rm m n, w 

1=0 m=-l 



where &; _m (i) = (— l) m 6™(t) for all to > (the over-bar denotes complex conjugation) and 6[j = 
(4-7r) -1 / 2 V t due to the normalization condition. Nematic states with fore-aft symmetry satisfy 
f(u) = /(— u), and for these I is restricted to the set of even integers. The macroscopic state of 
the suspension can be quantified by three variables - the structure tensor, S = (uu) — (5/3, the 
concomitant scalar structure factor S e = 9(S • S • S)/2] 1 / 3 and the mean polarity J = (u). We 



now specify the two inner products, {Y( n \f) = J*Y ; m (u) f(u,t) du, and mi\l2, m2\h, mz) = 
/ Y™ 1 (u) Y™ 2 (u) Y™ 3 (u) du and functions d 2n = [vr(4n+ l)(2n- 3)!!(2n- l)!!][2( 2n+2 )n!(n+ l)!]" 1 
and Co (l') = [(I' - - 3)!! 2 ][(/' + 2)(/'!!) 2 ]- 1 . 

Using these definitions with (4) we can write (2) as 

oo I' 

) 

2 v 15' 



l'=0 m '=-l 

and 

oo +21' 



V = -4kUY; E (JtT^' (u)b$ (6) 

i'=lm'=-2J' 1 f + j 

with {/ = 2NL 2 d. In writing (5) and (6) we have ignored constants linear in U and independent of 
u. The expressions are the same as those for non-magnetizable rods because the excluded volume 
potential is just dependent on geometrical symmetries. Parameters A' and £>' in (3) are proportional 
to the number density of the particles, and can be rewritten as A' = AU and B' = BU. Henceforth 
U, A and B are treated as three independent parameters. Combining (1), (4), (5) and (6) and 
using appropriate inner products we get the following evolution equation for the modes bf 1 , 

„ m oo +p 

_L. = _/(/ + i) hT _ £ £ {UEV + Um1 ( 7) 

p=0 q=—p 

where 

oo +/' 



-m=4^x; e ( 8 ) 

P=0m'=-I' 
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and oev depends on the nature of the excluded volume potential, 



ATT 00 + 1 ' 

VMS = ^Y1 E W'<^*, (9) 
l'=0 m'=-V 



oo +21' , 

do 



a =A,uY, E whw'*- do) 

l'=0 m'=-2V 

The function ^> = ^(Z, m, p, q, I', ml) is given by 

^(l,m,p,q,l',m') = —mm' (l,m\p, q\l' ,m') 

1 p(i + l)-m(m+l)] i , , 

2 ^ [/'(/' + 1) - m '(m' + l)]- 1 j ( ' + |P ' 91 ' + } 

2 [/'(/' + 1) — m (m' — I)] 1 

It is clear from equations (7)-(ll) that nematic branches corresponding to (A = 0, B > 0) and 
thus J = form a subset of possible stationary solutions to (7). It is also clear that (S = 0, J ^ 0) 
states are un-physical. 

A linear stability analysis of (7) about the isotropic state, / (u) = (47r) _1 is readily performed 
using b™ = (bf l )o+eb l rn +0(e 2 ), e«l being a suitable amplitude, and retaining terms through 0(e). 
The growth rates or eigenvalues, XT 1 , corresponding to the disturbance 5j m (u) can be obtained from 
the linearized equations. For the Maier-Saupe potential we get the following eigenvalues (for odd 
and even I respectively) (Af) M S = -1(1 + 1)(1 - <5j,i-4I7/3), and (Xf) M S = -1(1 + 1)(1 - U(l + 
^)<Jj2/5), indicating that there are two critical points on the 5 = isotropic branch. The first 
critical point satisfies (1 + B)U^ = 5. The critical eigenvalue is five fold degenerate with the 
associated destabilizing eigenvectors being linear combinations of Y™ ', m = —2,-1,0,1,2. The 
second critical point satisfies Jj\ = 3A" 1 and the critical eigenvalues that change sign at this point 
are three-fold degenerate and correspond to the eigenvectors Yj™, m = —1,0,1. In Figure (1) we 
plot these analytical predictions and compare them to numerically obtained solutions |l| for the 
case B = 1. We note that for fixed and finite B, as A — > oo, U\ — > 0. As A decreases from 
very large values, Jj\ < U£ initially and then, beyond a critical value of A, we get U b c > U%. 
For B = 1, the two critical points coincide for A = 1.2. Detailed numerical calculations show 
that for Uj! < U£, the branch is prolate, otherwise it is an oblate branch. For the Onsager 
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potential we find (for odd and even I respectively) (A™)o = — 1(1 + 1)(1 — AUS^i/3), and (A™)o = 
— + — C/(l + H)<5fe ) i/5 + C/7rc (Z)/2). Thus for odd I, as for the Maier-Saupe potential, there is 
one critical point on the S = line, [7* , which is the same as before. The destabilizing eigenvectors 
are the 3 independent components of Y 1 m (u). Let us denote the critical points for even I by U%(1) 
such that the critical eigenvectors at each point are the 21 + 1 independent components of YJ m (u). 
The first critical point occurs at U®(2) = (irc (2)/2 + B/5)" 1 and corresponds to the eigenvector 
set Y™(u). Higher order bifurcations occur at Uj!(l) = 2(irc (l)y 1 for I > 4 (k = 2,3, ..). 

We now concentrate on bifurcations of J > branches from the non-trivial J = nematic states 
for the specific case of a Maier-Saupe inter-molecular potential. As a point of departure to frame 
our discussion, we focus on the vicinity of the critical concentration given by = U b c and study 
the bifurcating branches as A and U are varied with B held fixed. 

Since the equations (1), (3), (7), (8) and (9) with .4 = exhibit rotational symmetry, we 
consider a base nematic state of the form (3) with coefficients {t>T)o rea l an d non-zero only if both 
I and m are even. Prom (1), (3) and (8) it is clear that the potential U and the parameter B can 
be combined into one dimensionless factor, W = U(l + B). Consider a base nematic state with 
director n = e 2 such that cos# = (u • n). Then the steady, uniaxial solution for this nematic is 
given by f(9) = exp (?>WS e cos 20/4) /P, where P is a normalizing constant. This yields 



plotted in Figure (2a). The solid lines are linearly stable branches. The oblate phase where the rods 
are oriented randomly in the (8 — nn) plane, is unstable to director fluctuations but stable if these 
are artificially suppressed - this is exemplified by the open circles which denote solutions obtained 
in integrating (1) in time in the subspace mentioned above 4 . Brownian dynamics simulations of the 
system for the Maier-Saupe potential 5 and B m = indicate that results using time integration for 
short times can yield an apparently stable oblate phase, thus mimicking for short times the effect 
of a pinned director. However long time integration of the stochastic system leads to the oblate 
branch being destabilized by symmetry breaking perturbations. We expect similar considerations 
to hold for B > 0. 

For later analysis we need an expression for the solution curve close to the critical point W = 5. 
An regular perturbation expansion in the small parameter, W = W — 5 indicates that along the 
nematic branches, we have the approximate relationship 



2S e + l 
3 




-i 



Jo z Jo 



Se(W) « -- 




+ o(w' A ) 



(12) 
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also plotted in Figure (2a) as the dash-dot line. We expect this to be accurate close to the critical 
point only. The structure factor for this nematic state has the form S D = — S e (W)S^ /3, with 
(s2x = Sy]j = —S^J/2). The eigenvalues obtained from (7) corresponding to the destabilizing 
eigenvectors, Y 2 m , are shown in Figure (2b). There are five eigenvalues that are zero at U%. The 
one corresponding to Y 2 ° (the structure parameter mode) has multiplicity of 1. The other four 
correspond to director fluctuations and occur as two pairs, one of which is identically zero. Since 
there are two independent ways to rotate a director on a sphere, we expect two neutral eigen- 
directions. 

We now impose small perturbations to the base state, b[ m , comprised only of even m modes while 
/ can be both even and odd. The equation for the growth of mode b'® with = 0, 2, 0, 1, 0) = 
and * (2) = (1, 0, 1, 0, 2, 0) = is: 

— --2^ (1 _ _ + __( &2 ) 0%) 



4ttU °° +P " 



-^EE Yj b'«(bf) o y{l,0,p,q,2,m')) (13) 

' p=0 q=-pm'=-2 

Close to criticality, the b'® mode dominates and so the p = 3 term in (13) can be ignored to leading 
order. Setting the growth rate to zero yields the following equation for A C (B, U) valid for small 

[! + _(!+ B)U(b 2 ) o * {2) ] = -^-(1 - 2tt(6S) * (1) ). (14) 

To obtain local information about the nature of the J > branches close to the critical point 
U£ = Uc, we expand all quantities in terms of a small parameter 5 that denotes the distance from the 
critical point measured along the (J = 0) nematic branches - to obtain (a) U = 5(l + £>) -1 (l + 5?7), 
(b) A c = 3(1 + + 5A c )/5 and (c) {b%) = 8{b\) w 5U' {d/ dU) Q {b\) = 5k m U with the slope 
Kn = — 7v / 5(10 v / 7r) _1 . Substituting these expressions in (14) yields at 0(6) 

A c = (27r(tf (1) + y [2) )k m - l)U = (15) 

Thus, close to the critical point as as we move along the prolate (with U locally decreasing), A c 
decreases as well. Similarly, as one moves along the oblate towards more higher values of U (U 
increases), A c increases. In short, critical points on the (J = 0, S e < 0) oblate state have A c > 1.2 
and on the (J = 0, S e > 0) prolate state satisfy A c < 1.2. 
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Our analysis yields insight about the behavior close to the critical point. Crucially, we find that 
it accords with numerical solutions far from the critical point obtained by Bhandar 1 for the specific 
case B = 1. Combining our local analytic results with these global numerical results, we obtain the 
bifurcation scenario illustrated in Figure 3. Let us recast the results in terms of the dependence of 
A c on the scalar structure parameter. For a fixed value of A, there are two critical points at which 
the J = branch becomes unstable to disturbances comprised of Yj° components. One of them 
is always on the S e = isotropic branch and the other is always on the (S e / 0, J = 0) nematic 
solution. When A < 1.2, the J > branches bifurcate at one point in the segment (S e = 0, 
U > 5/2) and at one point in the the prolate branch (J = , S e > 0. Even though the J = 
nematic prolate has a turning point at U « 2.245, the salient qualitative results of the local analysis 
holds even far from the critical point. 

Consider now the effects of an imposed external magnetic field H modeled by adding a term 
to the potential to (1) and (3) that is proportional to u • H. Such a field breaks the rotational 
degeneracy of the system inherent in (1). We anticipate that for a fixed values of U, A and B, 
the degree of order S as well as the extent of average polarization J change continuously with H. 
The transition from an isotropic to nematic state is replaced by a transition from a weakly aligned 
(paranematic) state to a strongly aligned state. Our results provide a mathematically convenient 
and physically relevant starting point to investigate these scenarios. 
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FIG. 1: A plot of the analytical value of A(U) for the Maier-Saupe potential at which the instability to 
Y™, m = — 1, 0, 1 modes arises on the isotropic J = S = branch. The circles are re-normalized computed 
results obtained from a numerical solution for B = 1 from Bhandar (2002) 1 . 
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(a) Bifurcation diagram when A=0 
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(b) Eigenvalues corresponding to critical eigenvectors emanating from U(1 +B)=5 



FIG. 2: (a) The equilibrium bifurcation diagram of the base nematic states with J = for „4 = 0. The 
prolate branch arising from U£ is unstable to structure factor fluctuations but regains stability beyond 
the turning point. The dash-dot line is the curve corresponding to the asymptotic expansion (12). (b) The 
eigenvalues corresponding to the destabilizing eigenvectors Y™ at U£ when ^4 = 0. The turning point is at 
W = U{1 + B) re 4.49 
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FIG. 3: Schematic sketch of the bifurcation scenario obtained by a combination of our local analytical 
results and global numerical results for B = 1. Region (A) corresponds to < U < U®, S e = J = and 
oo < A c < 1.2. As U increases, the critical value of A decreases, reaching 1.2 at U = U" = 5(1 + B)~ l . 
Region (B) corresponds to {/" < U < oo, S e = J = and 1.2 < A c < 0. Region (C) denotes bifurcation of 
(J > 0, S e > 0) nematic branches from the (J = 0, S e > 0) prolate curve. In this region, as one moves to 
S e — > 1, _4 C decreases from 1.2 to 0. Finally in region (D) along the oblate branch with (J = 0, S e < 0), we 
find A c increasing from 1.2 as S e decreases from to —1/2. 
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